#include "pointoperation.h"
#include <math.h>
#include <stdio.h>
#include <stdarg.h>

/**
 * 初始化点的坐标.
 *
 * \param point1 点
 */
void point_init(point point1)
{
	// 分配空间
	mpz_inits(point1->x, point1->y, NULL);

	// 初始化变量
	mpz_set_ui(point1->x, 0);
	mpz_set_ui(point1->y, 1);
}

/**
 * 初始化点的坐标.
 *
 * \param point1 点
 */
void point_inits(point_ptr point1, ...)
{
	// 预处理
	if (point1 == NULL)
	{
		return;
	}

	va_list valist;

	// 初始化 valist
	va_start(valist, point1);

	// 初始化第一个元素
	mpz_inits(point1->x, point1->y, NULL);
	mpz_set_ui(point1->x, 1);
	mpz_set_ui(point1->y, 1);

	// 依次获取元素，并初始化剩余元素
	point_ptr cur = va_arg(valist, point_ptr);

	while (cur != NULL)
	{
		mpz_inits(cur->x, cur->y, NULL);
		mpz_set_ui(cur->x, 1);
		mpz_set_ui(cur->y, 1);

		cur = va_arg(valist, point_ptr);
	}

	va_end(valist); // Clean up the va_list.
}

/**
 * 释放点的坐标所占的空间.
 *
 * \param point1 点
 */
void point_clear(point point1)
{
	mpz_clears(point1->x, point1->y, NULL);
}

/**
 * 释放点的元素所占的空间.
 *
 * \param point1 点
 * \param 可变参数 全是点+NULL
 */
void point_clears(point_ptr point1, ...)
{
	// 预处理
	if (point1 == NULL)
	{
		return;
	}

	va_list valist;

	// 初始化 valist
	va_start(valist, point1);

	// 初始化第一个元素
	mpz_clears(point1->x, point1->y, NULL);

	// 依次获取元素，并释放剩余元素的空间
	point_ptr cur = va_arg(valist, point_ptr);

	while (cur != NULL)
	{
		mpz_clears(cur->x, cur->y, NULL);

		cur = va_arg(valist, point_ptr);
	}
}

/**
 * 判断点是否值相等.
 *
 * \param point1 第一个点
 * \param point2 第二个点
 * \return x,y的值都相等，则返回0。有一个不相等，则返回-1
 */
int point_cmp(point point1, point point2)
{
	// x，y值都相等则返回0
	if (mpz_cmp(point1->x, point2->x) == 0 && mpz_cmp(point1->y, point2->y) == 0)
	{
		return 0;
	}

	return -1;
}

/**
 * 给点赋值.
 *
 * \param point1
 * \param point2
 */
void point_set(point point1, point point2)
{
	mpz_set(point1->x, point2->x);
	mpz_set(point1->y, point2->y);
}

/**
 * 初始化和赋值.
 *
 * \param point1
 * \param point2
 */
void point_init_set(point point1, point point2)
{
	// 分配空间
	mpz_inits(point1->x, point1->y, NULL);

	// 赋值
	mpz_set(point1->x, point2->x);
	mpz_set(point1->y, point2->y);
}

/**
 * 点加运算 ret = point1 + point2.
 *
 * \param point1 点P
 * \param point2 点Q
 * \param p 模数p
 * \param a 一次项系数
 * \return 点R
 */
void point_add(point ret, point point1, point point2, mpz_t a, mpz_t p)
{
	// 预处理 处理传过来的两个点、a和p
	mpz_t x1, y1, x2, y2, x3, y3, temp_a, temp_p;

	// 分配空间
	mpz_inits(x1, y1, x2, y2, x3, y3, temp_a, temp_p, NULL);

	// 变量赋值
	mpz_set(x1, point1->x);
	mpz_set(y1, point1->y);
	mpz_set(x2, point2->x);
	mpz_set(y2, point2->y);
	mpz_set_ui(x3, 0);
	mpz_set_ui(y3, 0);
	mpz_set(temp_a, a);
	mpz_set(temp_p, p);

	// 1.初始化lambda所需的分子、分母的变量
	mpz_t lambda, numerator, denominator, temp;
	// 分配空间
	mpz_inits(lambda, numerator, denominator, temp, NULL);
	// 赋值
	mpz_set_ui(lambda, 0);
	mpz_set_ui(numerator, 0);
	mpz_set_ui(denominator, 0);
	mpz_set_ui(temp, 0);

	// 2.计算lambda
	// lambda 有两种情况
	if (point_cmp(point1, point2) == 0)
	{
		// 第一个情况两个点相等 求导
		// (3x1^2+a)/(2y1) mod p
		// 2.1.1.计算分子 3x1^2+a
		// x1^2
		mpz_pow_ui(numerator, x1, 2);

		// 3x1^2
		mpz_set_ui(temp, 3);
		mpz_mul(numerator, temp, numerator);

		// 3x1^2+a
		mpz_add(numerator, numerator, temp_a);

		// 2.1.2.计算分母 2y1
		mpz_set_ui(temp, 2);
		mpz_mul(temp, temp, y1);

		// 2.1.3.分母求逆 (2y1)^(-1)
		inverse_by_p(denominator, temp, temp_p);

		// 2.1.4.计算出lambda (3x1^2+a)* (2y1)^(-1) mod p
		mpz_mul(lambda, numerator, denominator);
		mpz_mod(lambda, lambda, temp_p);
	}
	else
	{
		// 第二种情况两个点不相等
		// (y2-y1)/(x2-x1) mod p
		// 2.2.1.计算分子 (y2-y1) mod p
		mpz_sub(numerator, y2, y1);
		mpz_mod(numerator, numerator, temp_p);

		// 2.2.2.计算分母 (x2-x1) mod p
		mpz_sub(denominator, x2, x1);

		// 2.2.3.分母求逆 (x2-x1)^-1
		inverse_by_p(denominator, denominator, temp_p);

		// 2.2.4.计算出lambda (y2-y1)/(x2-x1) mod p
		mpz_mul(lambda, numerator, denominator);
		mpz_mod(lambda, lambda, temp_p);
	}

	// 3.通过计算出的lamda，计算x3，y3
	// 3.1.计算x3 x3=lambda^2-x1-x2 mod p
	// lambda^2
	mpz_pow_ui(x3, lambda, 2);
	// lambda^2-x1
	mpz_sub(x3, x3, x1);
	// lambda^2-x1-x2
	mpz_sub(x3, x3, x2);
	// lambda^2-x1-x2 mod p
	mpz_mod(x3, x3, temp_p);
	// 设置x3
	mpz_set(ret->x, x3);

	// 3.2.计算y3 y3=lambda(x1-x3)-y1 mod p
	// (x1-x3)  x3=ret->x
	mpz_sub(y3, x1, ret->x);
	// lambda(x1-x3)
	mpz_mul(y3, lambda, y3);
	// lambda(x1-x3)-y1
	mpz_sub(y3, y3, y1);
	// lambda(x1-x3)-y1 mod p
	mpz_mod(y3, y3, temp_p);
	// 设置y3
	mpz_set(ret->y, y3);

	// 4.善后处理
	// 释放空间
	mpz_clears(lambda, numerator, denominator, temp, NULL);
	mpz_clears(x1, y1, x2, y2, x3, y3, temp_a, temp_p, NULL);
}

/**
 * 点加的逆运算 ret = point1 - point2.
 *
 * \param point1 点P
 * \param point2 点Q
 * \param p 模数p
 * \param a 一次项系数
 * \return 点R
 */
void point_sub(point ret, point point1, point point2, mpz_t a, mpz_t p)
{
	// 1.先计算的 point2 的逆
	point point2_inv;
	point_init(point2_inv);
	// 设置 point2_inv->x = point2->x
	mpz_set(point2_inv->x, point2->x);

	mpz_t temp;
	mpz_init(temp);
	mpz_set_ui(temp, 0);

	// temp= -point2->y mod p;
	mpz_sub(temp, temp, point2->y);
	mpz_mod(temp, temp, p);
	// 设置 point2_inv->y = -point2->y
	mpz_set(point2_inv->y, temp);

	// 2.再计算 point1 + point2的逆
	point_add(ret, point1, point2_inv, a, p);

	// 3.释放空间
	mpz_clear(temp);
	point_clear(point2_inv);
}

/**
 * 多倍点运算。利用快速幂迭代法求解.
 *
 * \param ret 输出结果
 * \param point1 点
 * \param k 倍数
 * \param a 一次项系数
 * \param p 模数p
 */
void point_mul(point ret, point point1, mpz_t k, mpz_t a, mpz_t p)
{
	// 预处理
	// 拷贝传过来的输入
	point temp_point1;
	mpz_t temp_k, temp_a, temp_p;
	// 分配空间，并拷贝值
	point_init(temp_point1);
	mpz_set(temp_point1->x, point1->x);
	mpz_set(temp_point1->y, point1->y);

	mpz_inits(temp_k, temp_a, temp_p, NULL);
	mpz_set(temp_k, k);
	mpz_set(temp_a, a);
	mpz_set(temp_p, p);

	// 1.初始化
	// 用于输出 最开始是无穷远点 需要单独做判断
	point output_point;
	point_init(output_point);

	// 1表示是无穷远点
	int infinity_point = 1;

	// 常量0，1
	mpz_t const_zero, const_one, temp;
	mpz_inits(const_zero, const_one, temp, NULL);
	mpz_set_ui(const_zero, 0);
	mpz_set_ui(const_one, 1);
	mpz_set_ui(temp, 0);

	// 2.依次计算每一次的a^(2^i)
	while (mpz_cmp(temp_k, const_zero) != 0)
	{
		// 计算当前的2^i中的i是否为1，为1则与ret相乘
		// temp = b & 0x01;
		mpz_and(temp, temp_k, const_one);

		if (mpz_cmp(temp, const_one) == 0)
		{
			// 2.1. k & 0x01 = 1 ，则 ret = ret + point1

			// 最开始的时候，ret是无穷远点需要单独处理
			if (infinity_point == 1)
			{
				point_set(output_point, temp_point1);

				// ret 从无穷远点 变为 temp_point1
				// 因此 output_point的标志设置为 0
				infinity_point = 0;
			}
			else
			{
				point_add(output_point, output_point, temp_point1, temp_a, temp_p);
			}
		}
		// 2.2.更新temp_k
		// temp_k = temp_k / 2;
		mpz_tdiv_q_2exp(temp_k, temp_k, 1);

		// 2.3.更新 point1= 2 * point1
		point_add(temp_point1, temp_point1, temp_point1, temp_a, temp_p);
	}

	// 3.处理返回值
	point_set(ret, output_point);

	// 4.善后处理
	mpz_clears(const_zero, const_one, temp, NULL);

	point_clear(output_point);

	mpz_clears(temp_k, temp_a, temp_p, NULL);

	point_clear(temp_point1);
}

/**
 * 对 a^b mod p 求逆.
 *
 * \param ret 返回值
 * \param a 需要求逆元素的底数
 * \param b 需要求逆元素的指数
 * \param p 模数p
 * \return
 */
void inverse(mpz_t ret, mpz_t a, mpz_t b, mpz_t p)
{
	// 预处理
	// 用于拷贝传过来的a，b，p
	mpz_t temp_a, temp_b, temp_p;
	mpz_inits(temp_a, temp_b, temp_p, NULL);
	// 变量赋值
	mpz_set(temp_a, a);
	mpz_set(temp_b, b);
	mpz_set(temp_p, p);

	// a 应该小于 p
	mpz_mod(temp_a, temp_a, temp_p);
	// 返回值要做乘法计算，所以这里设置为1
	mpz_set_ui(ret, 1);

	// 1.初始化
	mpz_t zero, one, temp;
	// 分配空间
	mpz_inits(zero, one, temp, NULL);

	// 常数0
	mpz_set_ui(zero, 0);

	// 常数 1
	mpz_set_ui(one, 1);
	// 用于计算  b&0x01
	mpz_set_ui(temp, 0);

	// 依次计算每一次的2^i
	while (mpz_cmp(temp_b, zero) != 0)
	{
		// 计算当前的2^i中的i是否为1，为1则与ret相乘
		// temp = b & 0x01;
		mpz_and(temp, temp_b, one);
		if (mpz_cmp(temp, one) == 0)
		{
			// ret = (ret * temp_a) % p;
			mpz_mul(ret, ret, temp_a);
			mpz_mod(ret, ret, temp_p);
		}

		// 更新数据
		// 更新指数 (b>>1)
		// 右移动位运算
		// b= b / 2
		mpz_tdiv_q_2exp(temp_b, temp_b, 1);

		// 更新each
		// temp_a = (temp_a * temp_a) % p;
		mpz_mul(temp_a, temp_a, temp_a);
		mpz_mod(temp_a, temp_a, temp_p);
	}

	// 释放空间
	mpz_clears(zero, one, temp, NULL);
	mpz_clears(temp_a, temp_b, temp_p, NULL);

	// ret = ret % p;
	mpz_mod(ret, ret, temp_p);
}

/**
 * 费马小定理求逆。a^(p-2)=a^(-1) mod p.
 *
 * \param a 需要求逆得的元素
 * \param p 模数p
 * \return
 */
void inverse_by_p(mpz_t ret, mpz_t a, mpz_t p)
{
	mpz_t temp, two;
	mpz_inits(temp, two, NULL);
	mpz_set_ui(two, 2);

	// temp = p - 2;
	mpz_sub(temp, p, two);

	// ret= a^(p-2) % p
	inverse(ret, a, temp, p);

	// 释放空间
	mpz_clears(temp, two, NULL);
}

/**
 * 8.1.1 有限域上的椭圆曲线。计算椭圆曲线上的点.
 *
 * \param a 一次项系数
 * \param b 常数项
 * \param p 模数 mod为大于等于2的素数
 */
void cal_elliptic_curve_point(int a, int b, int p)
{
	// 输出提示
	printf("E%d(%d，%d)椭圆曲线点集\n", p, a, b);

	// 预处理
	int mod_judge_root = (p - 1) / 2;
	int mod_solve_root = (p + 1) / 4;
	// 1.x的取值范围是[0, p-1]
	for (int i = 0; i < p; i++)
	{
		// 2.计算出z=x^3+a^x+b
		int z = (i * i * i + a * i + b) % p;

		// 2.1.欧拉准则判别平方剩余是否有平方根
		int judge = pow(z, mod_judge_root);
		if (judge % p == 1)
		{
			// 2.1.1.有根
			// 解出两个点
			int positive_root = (int)pow(z, mod_solve_root) % p;
			int negative_root = (p - positive_root) % p;

			// 小的点在前面，大的点在后面
			int min = positive_root;
			int max = negative_root;
			if (min > max)
			{
				min = negative_root;
				max = positive_root;
			}

			printf("point(x,y) ：\t\t(%d，%d)\n", i, min);
			printf("point(x,y) ：\t\t(%d，%d)\n", i, max);
		}

		// 2.1.2.没有根 进行下一次循环
	}
}
